Load the data

library(tseries)
library(forecast)
library(TSstudio)

library(urca)
library(tidyquant)
library(magrittr)
library(MASS)
library(eFRED)

Rsxf <- fred("RSXFSN", all = FALSE)

2 Modification du data et détermination du type

Rsxf_ts = ts(Rsxf$RSXFSN, start = c(1992,01), frequency = 12)
ts_plot(Rsxf_ts, title = "Sales over time" )

On observe une trend avec une saisonalité qui se réperte et qui est additive car l’ampleur de la saisonnalité des fluctuations ne varient pas avec le niveau des séries chronologiques .

plot(decompose(Rsxf_ts))

On observe bien une trend et la sésonalité. On ne peut donc pas appliquer un ARMA(p,q) car on observe encore une non stationnarité et une sésonalité.

acf(Rsxf_ts, lag = 120)

pacf(Rsxf_ts, lag = 100 ) 

Notre Acf nous montre une grosse auto-corrélation donc dépendance des valeurs dans le passé ce qui montre une présence d’un AR On observe bien un ordre d’AR On observe aussi une saisonalité dans l’acf.

3 Filtrage saisonnalité

filter.outliers = function(timeseries)
{
library(foreach)
library(pracma)

## on décompose notre time series
stl_decomp = stl(timeseries, s.window = 'periodic', t.window = 13, robust = T)

## On remove la saisonnalité
TswithoutSais = seasadj(stl_decomp)

## apply a hampel filter
TsWS = hampel(TswithoutSais, k = 12)$y

## put back the seasonal component
return(TsWS)
}

Rsxf_Wseasonal <- filter.outliers(Rsxf_ts)

plot(Rsxf_Wseasonal)

On a enlevé la plus part de la saisonnalité.

#Autre méthode pour enlever la saisonnalité 
Rsxf_WithoutSesonal2 <- Rsxf_ts - decompose(Rsxf_ts)$seasonal
ts_plot(Rsxf_WithoutSesonal2)

plot(decompose(Rsxf_WithoutSesonal2))

Le même résultat

4 Check de notre Série

plot(decompose(Rsxf_Wseasonal)) 

On ne peut toujours pas utiliser notre signal car il reste encore la stationnarité même si on a enlevé la saisonnalité. Donc on ne peut toujours pas utilisé un modèle ARMA(p,q)

1 Etude stationnarité ADF

adf.test(Rsxf_Wseasonal)

## Warning in adf.test(Rsxf_Wseasonal): p-value greater than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Rsxf_Wseasonal
## Dickey-Fuller = 0.73695, Lag order = 7, p-value = 0.99
## alternative hypothesis: stationary

pvalue > O.O5 donc on ne rejette pas l’hypotese Ho qu’il y a présence de racine unitaire donc pi = 0 car (Pi = phi - 1) donc les coefficients sont significatif et Pi = 0 donc DS.

x <- ur.df(Rsxf_Wseasonal, type = "trend", lag = 2)
summary(x)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -58234  -5318    -37   5135  76858 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3501.96094 3483.86121   1.005    0.315    
## z.lag.1       -0.02160    0.02375  -0.910    0.364    
## tt            29.11833   23.56172   1.236    0.217    
## z.diff.lag1   -0.35715    0.05259  -6.791 4.58e-11 ***
## z.diff.lag2   -0.29702    0.05140  -5.779 1.62e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11220 on 362 degrees of freedom
## Multiple R-squared:  0.1729, Adjusted R-squared:  0.1638 
## F-statistic: 18.92 on 4 and 362 DF,  p-value: 3.835e-14
## 
## 
## Value of test-statistic is: -0.9096 5.0048 1.5297 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2  6.15  4.71  4.05
## phi3  8.34  6.30  5.36

#https://stats.stackexchange.com/questions/24072/interpreting-rs-ur-df-dickey-fuller-unit-root-test-results

On a pris lag 2 car il restait des * pour le lag = 1. Ici on voit que -0.9096 > -3.42 donc on ne rejette pas HO donc présence de racine unitaire donc \(\pi = 0\) donc la série n’est pas stationnaire. On ne rejette pas Ho quand Vtest > tau3 Ensuite 5.0048 > 4.71 ici on test Ho \(\beta1 = \beta 2 = \pi = 0\) on rejette donc Ho. On ne rejette pas Ho quand Vtest < phi2 Et enfin 1.5297 < 6.30 avec Ho \(\beta 2 = \pi = 0\) donc on ne rejette pas Ho. On ne rejette pas Ho quand Vtest < phi3

Donc la série n’est pas stationnaire avec $ = = 0 $ avec 5% d’erreur.

2 Degres d’intégration

D’après nos résultat, notre degrès d’intégration est de 2. Car le lag 1 et lag 2 sont significatifs pour le test ADF.

3 Phillips pheron

summary(ur.pp(Rsxf_Wseasonal))

## 
## ################################## 
## # Phillips-Perron Unit Root Test # 
## ################################## 
## 
## Test regression with intercept 
## 
## 
## Call:
## lm(formula = y ~ y.l1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -48064  -7039   -404   6371  92460 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.209e+03  2.018e+03   0.599    0.549    
## y.l1        1.000e+00  6.019e-03 166.131   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12320 on 367 degrees of freedom
## Multiple R-squared:  0.9869, Adjusted R-squared:  0.9868 
## F-statistic: 2.76e+04 on 1 and 367 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic, type: Z-alpha  is: 1.6289 
## 
##          aux. Z statistics
## Z-tau-mu           -0.1945

Ici pareil, on a une 1.6289 > -0.1945 donc on ne rejette pas H0 où \(\pi = 0\) et donc \(\phi =1\) présence de la racine unitaire et donc la série n’est pas stationnaire. On retrouve la même interprétation que Dickey et Fuller.

summary(ur.pp(diff(Rsxf_Wseasonal)))

## 
## ################################## 
## # Phillips-Perron Unit Root Test # 
## ################################## 
## 
## Test regression with intercept 
## 
## 
## Call:
## lm(formula = y ~ y.l1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -53527  -6369   -208   5529  85613 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1523.64961  619.45790   2.460   0.0144 *  
## y.l1          -0.28331    0.05019  -5.644 3.33e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11830 on 366 degrees of freedom
## Multiple R-squared:  0.08007,    Adjusted R-squared:  0.07756 
## F-statistic: 31.86 on 1 and 366 DF,  p-value: 3.328e-08
## 
## 
## Value of test-statistic, type: Z-alpha  is: -381.4512 
## 
##          aux. Z statistics
## Z-tau-mu            2.9123

Ici dans notre série filtrée de la trend et des saisonnalités, on observe donc bien que -381.4512 < 2.9123 donc on rejette H0 donc la série est stationnaire.

4 KPSS

z <- ur.kpss(Rsxf_Wseasonal, type = "tau")
summary(z)

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 5 lags. 
## 
## Value of test-statistic is: 0.5403 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216

Ici le test fait l’inverse donc 0.5403 > 0.146 donc on rejette l’hypothèse Ho donc la série n’est pas stationnaire. La série est donc bien non-stastionnaire.

5 Degres d’intégration avec KPSS

for (i in 1:7)
{
  x <- summary(ur.kpss(Rsxf_Wseasonal, type = "tau", use.lag = i))
  print(x)
}

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 1 lags. 
## 
## Value of test-statistic is: 1.4605 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216
## 
## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 2 lags. 
## 
## Value of test-statistic is: 1.0109 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216
## 
## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 3 lags. 
## 
## Value of test-statistic is: 0.7779 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216
## 
## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 4 lags. 
## 
## Value of test-statistic is: 0.6362 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216
## 
## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 5 lags. 
## 
## Value of test-statistic is: 0.5403 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216
## 
## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 6 lags. 
## 
## Value of test-statistic is: 0.4716 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216
## 
## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 7 lags. 
## 
## Value of test-statistic is: 0.4198 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216

Avec tous les test, on observe que \(\pi = 0\) et que \(\beta1\) est significatif. Cela correspond a une non-stationnarité processus DS Et que le degré d’intégration peut-être théoriquement de 2, or on voit que différentier avec un degré d’intégration de 1 rend notre série stationnaire.

6 on applique donc le filtre avec différention 1

Rsxf_Perfect <- diff(Rsxf_Wseasonal,lag = 1, differences = 1)
ts_plot(Rsxf_Perfect)

On ne voit plus de trend.

summary(ur.df(Rsxf_Perfect, type = "trend", lag = 1))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -57706  -5229    181   5063  76427 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  529.76076 1178.84973   0.449    0.653    
## z.lag.1       -1.68119    0.07973 -21.086  < 2e-16 ***
## tt             8.28801    5.54645   1.494    0.136    
## z.diff.lag     0.30842    0.04983   6.189 1.63e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11220 on 363 degrees of freedom
## Multiple R-squared:  0.6785, Adjusted R-squared:  0.6758 
## F-statistic: 255.3 on 3 and 363 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -21.0859 148.2183 222.3201 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2  6.15  4.71  4.05
## phi3  8.34  6.30  5.36

Ici -21.0859 << -3.42 donc rejet de Ho donc série stationnaire et donc \(\pi < 0\) Puis 148.2183 >> 4.71 donc non rejet de HO soit \(\beta 1 = \pi = 0\) Et enfin 222.3201 >> 6.30 donc non rejet de H0 soit \(\beta 1 = \beta 2 = \pi = 0\)

donc nous avons ici une série stationnaire avec un d = 1.

1 Most Revelant ARMA(p,q)

acf(Rsxf_Perfect, lag = 120)

On voit une petite autocorrélation avec les évènements passés donc peut-être présence d’un AR.

pacf(Rsxf_Perfect, lag = 120)

Avec la significativité différente de 0, on voit une corrélation différente de 0 pour potentiellement un AR(2)

On va tester le meilleur modèle par AIC et BIC :

best.order=c(0,0,0)
best.AIC=10000
for (i in 0:5) for (j in 0:5) 
{
  fit.AIC=AIC(arima(Rsxf_Perfect,order=c(i,0,j)))
  if(fit.AIC < best.AIC)
  {
    best.order=c(i,0,j)
    best.arma=arima(Rsxf_Perfect,order=best.order)
    best.AIC=fit.AIC
  }
}
 best.order 

## [1] 4 0 4

 best.AIC

## [1] 7894.76

best.order=c(0,0,0)
best.BIC= 10000
for (i in 0:5) for (j in 0:5)
{ 
  fit.BIC=BIC(arima(Rsxf_Perfect,order=c(i,0,j)))
  if(fit.BIC < best.BIC)
  {
    best.order=c(i,0,j)
    best.arma=arima(Rsxf_Perfect,order=best.order)
    best.BIC=fit.BIC
  }
}
best.order

## [1] 4 0 4

best.BIC

## [1] 7933.868

On obtient pour les 2 méthodes un ARMA(4,4).

Pour une alternative, on peut modéliser notre méthode avant le filtrage pour la rendre stationnaire avec un ARIMA(p,d,q), on va donc faire d’une pierre deux coups. Rendre la série stationnaire et lui trouver un bon model.

2 Test du model

Estimation <- arima(Rsxf_Perfect,order=c(4,0,4))

values_observed <- Rsxf_Perfect 

residus <- Estimation$residuals

values_estimated <-(values_observed - residus)

par(mfrow=c(2,1))

plot(values_estimated,type='l',ylab="values_estimated")
plot(Rsxf_Perfect, type = "l")

On voit une ressemblance même s’il le modèle estimé ne fit pas parfaitement.

Quality Checks : Residus

Quality_Checks = matrix(0,ncol=2,nrow=20)
for(i in 1:20)
{
  Quality_Checks[i,1] = Box.test(residus,lag=i,type="Ljung")$statistic
  Quality_Checks[i,2] = Box.test(residus,lag=i,type="Ljung")$p.value
}
plot(Quality_Checks[,1])

plot(Quality_Checks[,2])

On observe que la p-value > 0,05 les résidus sont donc indépendant sur cette intervalle. Cependant, Nous avons donc manqué quelque chose dans la modélisation de notre model. Peut-être le fait d’avoir fait une différence avec un degré d’uintégration à 1 au lieux de 2.

acf(residus)

pacf(residus)

On retrouve dans notre ACF, une petite autocorrélation mais ce modèle fit le mieux avec notre data.

qqnorm(residus)
qqline(residus)

On observe bien que les résidus suivent une distribution normale car la plupart des points se situent sur la ligne.

plot(density(residus))

On voit bien ici que les résidus trouvés forment bien une densité normale centré en 0.

QUALITY CHECKS Significativité des coefficients trouvés

\[RMSE = \sqrt{\frac{\sum_{i = 1}^{N} (Predicted_{i} - Actual_{i})^{2}}{N} }\]

RMSE <- function(x,y)
{
  # with x is the predicted i vector 
  # y the actual i vector
  rmse <- sqrt(mean((x-y)^2))
  return(rmse)
}

RMSE(values_estimated,values_observed)

## [1] 10395.05

1 Détermination degré d’intégration

summary(ur.df(jnj$close, type = "trend", lag = 2))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.3341  -2.0488  -0.0715   2.3043  14.5427 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.561670   0.560610   1.002 0.317193    
## z.lag.1     -0.021627   0.015865  -1.363 0.173828    
## tt           0.011900   0.007037   1.691 0.091858 .  
## z.diff.lag1 -0.094066   0.056652  -1.660 0.097861 .  
## z.diff.lag2 -0.208921   0.056455  -3.701 0.000255 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.163 on 304 degrees of freedom
## Multiple R-squared:  0.06368,    Adjusted R-squared:  0.05136 
## F-statistic: 5.169 on 4 and 304 DF,  p-value: 0.0004839
## 
## 
## Value of test-statistic is: -1.3632 3.3917 1.5724 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2  6.15  4.71  4.05
## phi3  8.34  6.30  5.36

On observe bien que -1,3 > -3,42 donc on ne rejette pas H0 soit \(\pi = 0\) donc série non stationnaire. J’ai fais la l’ordre 2 car on a une p-value tres significative 3* avec une degré 2. On suppose donc un degré d’intégration de 2.

summary(ur.df(diff(jnj$close), type = "trend", lag = 1))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.6654  -2.0530  -0.0332   2.2828  14.3702 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.168653   0.477828   0.353   0.7244    
## z.lag.1     -1.327957   0.082388 -16.118   <2e-16 ***
## tt           0.003020   0.002666   1.133   0.2582    
## z.diff.lag   0.220373   0.055905   3.942   0.0001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.169 on 305 degrees of freedom
## Multiple R-squared:  0.5663, Adjusted R-squared:  0.562 
## F-statistic: 132.7 on 3 and 305 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -16.1183 86.6028 129.9037 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2  6.15  4.71  4.05
## phi3  8.34  6.30  5.36

Or ici j’intègre avec un degré de 1 et cela me renvoie a une série stationnaire car -16 < -3 donc on rejette H0 => donc série stationnaire.

Au final le degré d’intégration est de 1.

2 Order ARIMA(p,d,q)

On connait d = 1 donc on peut faire la méthode du AIC et BIC pour déterminer l’ordre du ARIMA.

best.order=c(0,1,0)
best.AIC=10000
for (i in 0:5) for (j in 0:5) 
{
  fit.AIC=AIC(arima(jnj$close,order=c(i,1,j)))
  if(fit.AIC < best.AIC)
  {
    best.order=c(i,1,j)
    best.arma=arima(jnj$close,order=best.order)
    best.AIC=fit.AIC
  }
}
 best.order 

## [1] 4 1 5

 best.AIC

## [1] 1772.018

best.order=c(0,1,0)
best.BIC= 10000
for (i in 0:5) for (j in 0:5)
{ 
  fit.BIC=BIC(arima(jnj$close,order=c(i,1,j)))
  if(fit.BIC < best.BIC)
  {
    best.order=c(i,1,j)
    best.arma=arima(jnj$close,order=best.order)
    best.BIC=fit.BIC
  }
}
best.order

## [1] 0 1 2

best.BIC

## [1] 1789.379

On se retrouve avec 2 modeles possible soit : AIC ARIMA(5,1,5) BIC ARIMA(0,1,2)

Cependant, avec notre PACF nous avions vu qu’il n’y avait pas d’ordre pour les AR donc le BIC semble être le plus représentatif.

auto.arima(jnj$close)

## Series: jnj$close 
## ARIMA(2,1,2) with drift 
## 
## Coefficients:
##           ar1     ar2     ma1      ma2   drift
##       -0.2286  0.4235  0.0977  -0.6762  0.4724
## s.e.   0.1591  0.1557  0.1329   0.1322  0.1225
## 
## sigma^2 = 17.05:  log likelihood = -879.94
## AIC=1771.87   AICc=1772.15   BIC=1794.31

Avec une auutre méthode, on retrouve le modèle ARIMA(0,1,2) comme étant le plus significatif.

3 Fit ARIMA estimated valued

Estimation <- arima(jnj$close,order=c(0,1,2))

values_observed <- jnj$close

residus <- Estimation$residuals

values_estimated <-(values_observed - residus)

par(mfrow=c(2,1))

plot(values_estimated,type='l',ylab="values_estimated")
plot(jnj$close, type = "l")

On observe une grosse ressemblence.

4 QUality checks on Residuals

Quality_Checks = matrix(0,ncol=2,nrow=20)
for(i in 1:10)
{
  Quality_Checks[i,1] = Box.test(residus,lag=i,type="Ljung")$statistic
  Quality_Checks[i,2] = Box.test(residus,lag=i,type="Ljung")$p.value
}
plot(Quality_Checks[,1])

plot(Quality_Checks[,2])

On observe que la p-value > 0,05 les résidus sont donc indépendant sur cette intervalle. Cependant, on a un drop de pvalue donc nous avons donc manqué quelque chose dans la modélisation de notre model. Peut-être le fait d’avoir fait une différence avec un degré d’uintégration à 1 au lieux de 2.

acf(residus)

#pacf(residus)

On voit dans l’ACF de nos résidus sont indépendants les uns des autres.

qqnorm(residus)
qqline(residus)

On observe bien que les résidus suivent une distribution normale car la plupart des points se situent sur la ligne.

plot(density(residus))

On voit bien ici que les résidus trouvés forment bien une densité normale centré en 0.

5 Forecast

predic <-forecast(values_estimated,level=c(95),h=10)

plot(predic, xlim=c(250,320))
prediction<- predict(predic,jnj$close,se = TRUE,interval="prediction")
lines(prediction$lower,lty="dashed",col="red")
lines(prediction$upper,lty="dashed",col="red")

CALULER CONFIDENCE DU FORCAST

1 Zivot and Andrews

Le test de Zivot et Andrews améliore le test de stationnarité d’une série temporelle. En effet, cette méthode dait un test de racine unitaire avec une rupture structurelle. Cette rupture est estimé.

Soit notre Hypothèse H0 : \(\pi = 0\) donc série non stationnaire H1 : \(\pi < 0\) avec soit : - 0 break - 1 break -> constante -> Trend -> constante + Trend

2 Generer 3 marches aléatoires

set.seed(123)
ut <- rnorm(500,0,1)

t <- 1:500
# Random walk
y <- cumsum(ut)

# Random walk with break 
ybreak <- c()
for (i in 1:500)
{
  ifelse( i< 250, ybreak[i] <- cumsum(ut[i]),ybreak[i] <- 20 + cumsum(ut[i]))
}

#random walk with break and trend
ytrend <- c()
for (i in 1:500)
{
  ifelse( i< 250, ytrend[i] <- cumsum(ut[i]) + 0.05*t[i] ,ytrend[i] <- 20 + cumsum(ut[i]) + 0.05*t[i])
}

plot(ytrend, type = "l")
lines(ybreak, type = "l", col = "blue")
lines(y, type = "l", col = "red")

#library(strucchange)
#library(zoo)

#nbreak <- breakpoints(ybreak~t, h = 20) 

3 Test Zivot and Andrews

Dans un premier temps on va regarder avec un test adf avec comme condition On ne rejette pas Ho quand Vtest > tau3 donc série non stationnaire.

#summary(ur.df(y, type = "trend", lag = 1))
#summary(ur.df(ybreak, type = "trend", lag = 1))
#summary(ur.df(ytrend, type = "trend", lag = 1))

# J'évite de le print mais je vous met une interpretation en dessous.

Pour y la marche aléatoire normale : -2.41 > -3.42 donc la série n’est pas stationnaire De même pôur la marche aléatoire avec un break (constante) : -2.6139 > -3.42 Et pour la marche aléatoire avec break et trend -2.6139 > -3.42

Donc d’après notre test ADF aucune des 3 n’est stationnaires.

Maintennant on applique le test de Zivot et Andrews.

summary(ur.za(y, model = "both", lag = 1))

## 
## ################################ 
## # Zivot-Andrews Unit Root Test # 
## ################################ 
## 
## 
## Call:
## lm(formula = testmat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.95301 -0.59372 -0.02188  0.61290  2.99265 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.2298772  0.1288764   1.784  0.07509 .  
## y.l1         0.9486190  0.0140338  67.595  < 2e-16 ***
## trend       -0.0012935  0.0008162  -1.585  0.11366    
## y.dl1       -0.0347155  0.0448744  -0.774  0.43953    
## du           0.4979822  0.1917264   2.597  0.00968 ** 
## dt           0.0023018  0.0013215   1.742  0.08217 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.963 on 492 degrees of freedom
##   (2 observations effacées parce que manquantes)
## Multiple R-squared:  0.9613, Adjusted R-squared:  0.9609 
## F-statistic:  2445 on 5 and 492 DF,  p-value: < 2.2e-16
## 
## 
## Teststatistic: -3.6612 
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82 
## 
## Potential break point at position: 263

summary(ur.za(ybreak, model = "both", lag = 1))

## 
## ################################ 
## # Zivot-Andrews Unit Root Test # 
## ################################ 
## 
## 
## Call:
## lm(formula = testmat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7389 -0.6166 -0.0155  0.6387  3.2192 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.0736376  0.1261368   0.584    0.560    
## y.l1        -0.0444094  0.0363614  -1.221    0.223    
## trend       -0.0006451  0.0008713  -0.740    0.459    
## y.dl1        0.0092039  0.0287074   0.321    0.749    
## du          21.0483874  0.7376339  28.535   <2e-16 ***
## dt           0.0006860  0.0012232   0.561    0.575    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9762 on 492 degrees of freedom
##   (2 observations effacées parce que manquantes)
## Multiple R-squared:  0.9908, Adjusted R-squared:  0.9907 
## F-statistic: 1.054e+04 on 5 and 492 DF,  p-value: < 2.2e-16
## 
## 
## Teststatistic: -28.723 
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82 
## 
## Potential break point at position: 249

summary(ur.za(ytrend, model = "both", lag = 1))

## 
## ################################ 
## # Zivot-Andrews Unit Root Test # 
## ################################ 
## 
## 
## Call:
## lm(formula = testmat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7389 -0.6166 -0.0155  0.6387  3.2192 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.070957   0.126134   0.563    0.574    
## y.l1        -0.044409   0.036361  -1.221    0.223    
## trend        0.051575   0.002001  25.780   <2e-16 ***
## y.dl1        0.009204   0.028707   0.321    0.749    
## du          21.048387   0.737634  28.535   <2e-16 ***
## dt           0.000686   0.001223   0.561    0.575    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9762 on 492 degrees of freedom
##   (2 observations effacées parce que manquantes)
## Multiple R-squared:  0.9966, Adjusted R-squared:  0.9966 
## F-statistic: 2.899e+04 on 5 and 492 DF,  p-value: < 2.2e-16
## 
## 
## Teststatistic: -28.723 
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82 
## 
## Potential break point at position: 249

Pour la marche aléatoire simple : -3.6612 > -5.08 donc non rejet de H0 donc \(\pi = 0\) donc non stationaire avec potentiellement un break à l’indice 263.

Cependant pour la marche aléatoire avec break (en constante) : -28.723 < -5.08 donc on rejette H0 donc la série est stationnaire avec un break à l’indice 249, pour rappelle j’ai mis mon break à l’indice 250. Le test est donc assez précis.

Et enfin pour la marche aléatoire avec break et trend, on a -28.723 < -5.08 donc de même série stationnaire mais avec un break à l’indice 249.

on voit donc l’importance de faire des tests avec ruptures structurelles.

4 ZA in US retail sales

on replot notre jeu de donnée sans séonalité.

ts_plot(Rsxf_Wseasonal)

On observe potentiellement un break au alentour des années 2007.

#summary(ur.df(Rsxf_Wseasonal, type = "trend", lag = 1))
summary(ur.za(Rsxf_Wseasonal, model = "both", lag = 1))

## 
## ################################ 
## # Zivot-Andrews Unit Root Test # 
## ################################ 
## 
## 
## Call:
## lm(formula = testmat)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -57499  -6186    314   5521  54650 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.251e+04  6.053e+03   7.023 1.08e-11 ***
## y.l1         7.176e-01  4.024e-02  17.833  < 2e-16 ***
## trend        2.515e+02  3.605e+01   6.977 1.44e-11 ***
## y.dl1       -1.607e-01  5.096e-02  -3.153  0.00175 ** 
## du           3.749e+04  5.807e+03   6.457 3.45e-10 ***
## dt          -2.886e+02  4.486e+02  -0.643  0.52048    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10940 on 362 degrees of freedom
##   (2 observations effacées parce que manquantes)
## Multiple R-squared:  0.9897, Adjusted R-squared:  0.9896 
## F-statistic:  6972 on 5 and 362 DF,  p-value: < 2.2e-16
## 
## 
## Teststatistic: -7.0187 
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82 
## 
## Potential break point at position: 350

Je rappelle qu’avec DF on trouve : -2.3639 > -3.42 donc séries non stationnaire Alors qu’avec ZA on trouve : -7.0187 < -5.08 donc série stationnaire avec un break à l’indice 350 Donc on peut se poser la question si ce break, chute en 2007, à l’indice 350 ne faussait pas tous nos tests.

Load data

Gdp <- fred("GDP", all = FALSE)
Gdp <- na.omit(Gdp, "GDP")

J’enlève toute valeure null.

Modification en TS

Gdp_ts = ts(Gdp$GDP, start = c(1947,01), frequency = 4)
ts_plot(Gdp_ts, title = "Gross Domestic Product over time" )

On observe un belle trend avec pas de saisonalité apparente.

Tests sur notre Time serie

#plot(decompose(Gdp_ts))

acf(Gdp_ts, lag = 100)

pacf(Gdp_ts)

ACF -> Grosse autocorrélation, peut-être présence d’un AR. PACF -> 1 seul corrélation à 1

summary(ur.df(Gdp_ts, type = "trend", lag = 1))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2190.79   -10.84    -1.27    19.13  1231.84 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.007254  26.986428   0.000   0.9998  
## z.lag.1      0.009917   0.004156   2.386   0.0176 *
## tt           0.184980   0.323210   0.572   0.5675  
## z.diff.lag  -0.128714   0.058249  -2.210   0.0279 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 170.4 on 297 degrees of freedom
## Multiple R-squared:  0.1706, Adjusted R-squared:  0.1623 
## F-statistic: 20.37 on 3 and 297 DF,  p-value: 4.955e-12
## 
## 
## Value of test-statistic is: 2.3865 38.4439 30.0376 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2  6.15  4.71  4.05
## phi3  8.34  6.30  5.36

On observe avec le test de Dickey Fuller que 2.3865 > -3.42 donc non rejet de H0 -> \(\pi = 0\) donc série non stationnaire. Par rapport au drift et trend notre T-stat sont supérieures aux valeurs 5% donc on rejette leur H0 donc rejet de nullité.

Mais nous aller voir si ce n’est pas un problème de break

summary(ur.za(Gdp_ts, model = "both", lag = 1))

## 
## ################################ 
## # Zivot-Andrews Unit Root Test # 
## ################################ 
## 
## 
## Call:
## lm(formula = testmat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2078.14   -21.93     0.01    30.15   380.78 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -34.196952  23.047319  -1.484 0.138937    
## y.l1          0.995171   0.003762 264.545  < 2e-16 ***
## trend         0.969200   0.282956   3.425 0.000701 ***
## y.dl1        -0.178682   0.051268  -3.485 0.000566 ***
## du          909.832281 106.208351   8.566 6.03e-16 ***
## dt          -61.533073  19.202554  -3.204 0.001502 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 141.9 on 295 degrees of freedom
##   (2 observations effacées parce que manquantes)
## Multiple R-squared:  0.9996, Adjusted R-squared:  0.9996 
## F-statistic: 1.439e+05 on 5 and 295 DF,  p-value: < 2.2e-16
## 
## 
## Teststatistic: -1.2838 
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82 
## 
## Potential break point at position: 294

On a -1.2838 > -5.08 donc non rejet de H0 qui comme adf \(\pi = 0\) donc série non sationnaire mais potentielle break à l’indice 294.

Conclusion : On sait donc que notre série est bien non sattionnaire avec d = 2

Différentiation et estimation du model ARIMA

Différention avec d = 1

Gdp_perfect <- diff(Gdp_ts, differences = 2)
ts_plot(Gdp_perfect)

Test de stationnarité

summary(ur.df(Gdp_perfect, type = "trend", lag = 1))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2055.74   -18.99     1.60    14.27  2067.12 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.37458   22.33321  -0.151    0.880    
## z.lag.1     -2.19791    0.09628 -22.828  < 2e-16 ***
## tt           0.04451    0.12842   0.347    0.729    
## z.diff.lag   0.37815    0.05392   7.013 1.59e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 191.6 on 295 degrees of freedom
## Multiple R-squared:  0.8264, Adjusted R-squared:  0.8246 
## F-statistic:   468 on 3 and 295 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -22.8284 173.7121 260.5682 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2  6.15  4.71  4.05
## phi3  8.34  6.30  5.36

#summary(ur.pp(Gdp_perfect))
#summary(ur.za(Gdp_perfect, model = "both", lag = 1))

avec notre test on a -22.8284 < -3.42 donc rejet de H0 donc série stationnaire.

Estimation du model ARIMA avec d = 2 et ARMA

acf(Gdp_perfect)

pacf(Gdp_perfect)

ACF -> Très petit autocorrélation mais existante MA(1) PACF -> petite corrélation AR(1)

auto.arima(Gdp_ts)

## Series: Gdp_ts 
## ARIMA(0,2,3) 
## 
## Coefficients:
##           ma1     ma2      ma3
##       -1.0690  0.2205  -0.0932
## s.e.   0.0579  0.0859   0.0577
## 
## sigma^2 = 30054:  log likelihood = -1978.49
## AIC=3964.97   AICc=3965.11   BIC=3979.8

auto.arima(Gdp_perfect)

## Series: Gdp_perfect 
## ARIMA(0,0,2) with non-zero mean 
## 
## Coefficients:
##           ma1     ma2    mean
##       -1.0700  0.0990  0.9300
## s.e.   0.0533  0.0575  0.3959
## 
## sigma^2 = 29847:  log likelihood = -1977.81
## AIC=3963.61   AICc=3963.75   BIC=3978.44

ARIMA(p,2,q)

best.order=c(0,2,0)
best.AIC=10000
for (i in 0:5) for (j in 0:5) 
{
  fit.AIC=AIC(arima(Gdp_ts,order=c(i,2,j), method = "ML"))
  if(fit.AIC < best.AIC)
  {
    best.order=c(i,2,j)
    best.arma=arima(Gdp_ts,order=best.order, method = "ML")
    best.AIC=fit.AIC
  }
}
 best.order 

## [1] 1 2 1

 best.AIC

## [1] 3964.754

best.order=c(0,2,0)
best.BIC= 10000
for (i in 0:5) for (j in 0:5)
{ 
  fit.BIC=BIC(arima(Gdp_ts,order=c(i,2,j), method = "ML"))
  if(fit.BIC < best.BIC)
  {
    best.order=c(i,2,j)
    best.arma=arima(Gdp_ts,order=best.order, method = "ML")
    best.BIC=fit.BIC
  }
}
best.order

## [1] 0 2 1

best.BIC

## [1] 3975.57

Avec le AIC, on obtient un ARIMA(1,2,1) et pour le BIC un ARIMA(0,2,1)

Determination ARMA(p,q) avec notre série différentié

best.order=c(0,0,0)
best.AIC=10000
for (i in 0:5) for (j in 0:5) 
{
  fit.AIC=AIC(arima(Gdp_perfect,order=c(i,0,j), method = "ML"))
  if(fit.AIC < best.AIC)
  {
    best.order=c(i,0,j)
    best.arma=arima(Gdp_perfect,order=best.order, method = "ML")
    best.AIC=fit.AIC
  }
}
 best.order 

## [1] 0 0 3

 best.AIC

## [1] 3961.523

best.order=c(0,0,0)
best.BIC= 10000
for (i in 0:5) for (j in 0:5)
{ 
  fit.BIC=BIC(arima(Gdp_perfect,order=c(i,0,j), method = "ML"))
  if(fit.BIC < best.BIC)
  {
    best.order=c(i,0,j)
    best.arma=arima(Gdp_perfect,order=best.order, method = "ML")
    best.BIC=fit.BIC
  }
}
best.order

## [1] 0 0 1

best.BIC

## [1] 3975.339

Avec le AIC, on obtient un ARIMA(0,0,3) et pour le BIC un ARIMA(0,0,1)

En conclusion, j’utilise le ARIMA(1,2,1) car il est correct par rapport à notre ACF et PACF Et l’ARMA(0,0,2)

Valeurs estimées avec nos 2 modèles ARIMA(1,2,1) et ARMA(0,0,2)

ARIMA(1,2,1)

Estimation_Gdp <- arima(Gdp_ts,order=c(1,2,1), method = "ML")

values_observed_Gdp <- Gdp_ts

residus_Gdp <- Estimation_Gdp$residuals

values_estimated_Gdp <-(values_observed_Gdp - residus_Gdp)

par(mfrow=c(2,1))

plot(values_estimated_Gdp,type='l',ylab="values_estimated")
plot(Gdp_ts, type = "l")

ARMA(0,0,2)

Estimation_Gdp2 <- arima(Gdp_perfect,order=c(0,0,2), method = "ML")

values_observed_Gdp2 <- Gdp_perfect

residus_Gdp2 <- Estimation_Gdp2$residuals

values_estimated_Gdp2 <-(values_observed_Gdp2 - residus_Gdp2)

par(mfrow=c(2,1))

plot(values_estimated_Gdp2,type='l',ylab="values_estimated")
plot(Gdp_perfect, type = "l")

Quality checks

ARIMA(1,2,1)

Quality_Checks = matrix(0,ncol=2,nrow=20)
for(i in 1:10)
{
  Quality_Checks[i,1] = Box.test(residus_Gdp,lag=i,type="Ljung")$statistic
  Quality_Checks[i,2] = Box.test(residus_Gdp,lag=i,type="Ljung")$p.value
}
plot(Quality_Checks[,1])

plot(Quality_Checks[,2])

On observe que la p-value > 0,05 jusqu’à un lag < 10 les résidus sont donc indépendant sur cette intervalle. Nous avons donc manqué quelque chose dans la modélisation de notre model.

acf(residus_Gdp)

#pacf(residus_Gdp)

On voit dans l’ACF de nos résidus sont indépendants les uns des autres.

qqnorm(residus_Gdp)
qqline(residus_Gdp)

On observe bien que les résidus suivent une distribution normale car la plupart des points se situent sur la ligne.

plot(density(residus_Gdp))

On voit bien ici que les résidus trouvés forment bien une densité normale centré en 0.

RMSE(values_estimated_Gdp,values_observed_Gdp)

## [1] 172.4497

ARMA(0,0,2)

Quality_Checks = matrix(0,ncol=2,nrow=20)
for(i in 1:10)
{
  Quality_Checks[i,1] = Box.test(residus_Gdp2,lag=i,type="Ljung")$statistic
  Quality_Checks[i,2] = Box.test(residus_Gdp2,lag=i,type="Ljung")$p.value
}
plot(Quality_Checks[,1])

plot(Quality_Checks[,2])

On observe que la p-value > 0,05 les résidus sont donc indépendant sur cette intervalle.

acf(residus_Gdp2)

#pacf(residus_Gdp2)

On voit dans l’ACF de nos résidus sont indépendants les uns des autres.

qqnorm(residus_Gdp2)
qqline(residus_Gdp2)

On observe bien que les résidus suivent une distribution normale car la plupart des points se situent sur la ligne.

plot(density(residus_Gdp2))

On voit bien ici que les résidus trouvés forment bien une densité normale centré en 0.

RMSE(values_estimated_Gdp2,values_observed_Gdp2)

## [1] 171.9032

On remarque que les résidus des 2 modèles sont très similaires en terme de leur forme et la même RMSE

Forecasting

ARIMA(1,2,1)

predic_Gdp <-forecast(values_estimated_Gdp,level=c(95),h=10)

plot(predic_Gdp, xlim=c(2000,2025))
prediction_Gdp<- predict(predic_Gdp,Gdp_ts,se.fit = TRUE,interval="prediction")
lines(prediction_Gdp$lower,lty="dashed",col="red")
lines(prediction_Gdp$upper,lty="dashed",col="red")

ARMA(0,0,2)

predic_Gdp2 <-forecast(Estimation_Gdp2,level=c(95),h=5)

plot(predic_Gdp2,xlim=c(2015,2024))
prediction_Gdp2<- predict(predic_Gdp2,Gdp_perfect,se.fit = TRUE, interval="prediction")
lines(prediction_Gdp2$lower,lty="dashed",col="red")
lines(prediction_Gdp2$upper,lty="dashed",col="red")

Voici les différents forecasting. On remarque que le forecasting du model ARIMA(1,2,1) semble meilleur que celui de ARMA(0,0,2).